Tensor Renormalization Group: Local Magnetizations, Correlation Functions, and 
Phase Diagrams of Systems with Quenched Randomness 
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The tensor renormalization-group method, developed by Levin and Nave, brings systematic im- 
provability to the position-space renormalization-group method and yields essentially exact results 
for phase diagrams and entire thermodynamic functions. The method, previously used on systems 
with no quenched randomness, is extended in this study to systems with quenched randomness. Lo- 
cal magnetizations and correlation functions as a function of spin separation are calculated as tensor 
products subject to renormalization-group transformation. Phase diagrams are extracted from the 
long-distance behavior of the correlation functions. The approach is illustrated with the quenched 
bond-diluted Ising model on the triangular lattice. An accurate phase diagram is obtained in tem- 
perature and bond-dilution probability, for the entire temperature range down to the percolation 
threshold at zero temperature. 
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I. INTRODUCTION 

The tensor renormalization-group (TRG) method de- 
veloped by Levin and Nave [l[ is a highly useful up- 
date of the traditional position-space renormalization- 
group approaches. While these founding approaches 
relied on uncontrolled approximations that were often 
system-specific 0-0], the TRG is general in scope — it 
works on any classical two-dimensional lattice Hamil- 
tonian with local interactions — and its accuracy can be 
systematically improved to converge on the exact ther- 
modynamic results. Along with these advantages, the 
method fits within the conceptual framework of tradi- 
tional renormalization-group theory: it is a mapping be- 
tween Hamiltonians on the original and coarse-grained 
lattices, and phase transition behavior can be extracted 
from flows of the Hamiltonians as the transformation is 
iterated 

The initial TRG study demonstrated the power of the 
approach in the context of the triangular-lattice Ising 
model Since then it has proven a versatile tool for 
a variety of classical systems, including the frustrated 
Ising model on a Shastry-Sutherland lattice [9\ , relevant 
to magnetization plateaus in rare-earth tetraborides, and 
the zero-hopping limit of a model for ultra-cold bosonic 
polar molecules on a hexagonal optical lattice • More- 
over, the ideas behind the TRG method have become the 
kernel for developments in two-dimensional quantum sys- 
tems pdl - [l6| . most notably tensor-entanglement renor- 
malization group for studying symmetry breaking and 
topological phase transitions [11], and accurate methods 
to calculate ground-state expectation values [l2l - fl4j . Be- 
yond the precision of the method, a key factor spurring 
the growth of tensor RG applications in both classical 
and quantum cases is computational efficiency: the CPU 
cost of carrying out TRG scales linearly with lattice 



size 



14] 



Given these promising characteristics, TRG is a 
natural candidate for tackling models with quenched 
randomness — a field where extracting accurate phase di- 
agram information is a significant challenge. The current 
study presents the first example of TRG applied to such 
a system with frozen disorder, namely the percolative 
system of the bond-diluted triangular-lattice Ising fer- 
romagnet, yielding, as seen in Fig.l, a highly accurate 
global phase diagram, down to zero temperature, where 
it connects with the percolation transition. 

Our paper is organized as follows: Sec. [TTJ develops 
the TRG method for a general quenched random system. 
Sec. lIIII illustrates this tensor network mapping in partic- 
ular for the bond-diluted model and shows how to extract 
physical observables such as spin-spin correlation func- 
tions. Sec. II VI uses this method, together with finite-size 
scaling relations for the correlation functions, to derive 
our main result: The phase diagram in terms of tempera- 
ture vs. bond dilution probability. Close agreement with 
the known critical temperature curve fi~7| is achieved even 
at a relatively low order of the TRG approximation (i.e., 
a small cutoff parameter). Our work opens up future 
possibilities for the extensive use of TRG in quenched 
disordered systems, as argued in the concluding remarks 
of Sec. El 



II. TRG METHOD FOR QUENCHED RANDOM 
SYSTEMS 



A. The Tensor Network 

As in earlier studies [H, we focus here on classi- 
cal Hamiltonians associated with hexagonal-lattice ten- 
sor networks, though the method that we develop for 



2 



4.0 



3.0 

i — i 

01 
i_ 

ro 2.0 

Q_ 



1.0 



Exact 

• D=12 



Para 



Ferro 



n=0 



n=l 



n=2 



0.2 0.4 0.6 0.8 

Bond dilution probability p 



1.0 



FIG. 1: (Color online) The phase diagram of the bond-diluted 
Ising model on a triangular lattice, showing the transition 
temperature as a function of the bond dilution probability p. 
The ferromagnetic (Ferro) and paramagnetic (Para) phases 
are marked. The phase boundary line between these two 
phases connects, at zero temperature, with the percolation 
transition on the triangular lattice. Filled circles are our re- 
sults using the TRG method with D = 12 together with finite 
size scaling, as described in Sec. IIVI The red dotted line is 
the result of Georges et al. [13], which is exact on the scale of 
the figure. 



quenched random systems is readily generalized to other 
geometries like the square and kagome lattices [![. We 
consider a general Hamiltonian that involves local inter- 
actions expressed in terms of bond degrees of freedom, 
such that each bond has d possible states and the parti- 
tion function of the system has the form 



ii,...,iK=l 



■T, 



(1) 



where, for each of the N sites in the hexagonal lattice, 
a real- valued tensor Ti m i„i is a Boltzmann weight de- 
pending on the configuration of the three bonds meeting 
at the site. The bond degrees of freedom correspond to 
each tensor index running from 1 to d. These bond in- 
dices are labeled i\ through ix, for the total of K = 3N/2 
bonds in the lattice. Although the tensor can have as 
many as d 3 distinct non-zero elements, in practice some 
bond configurations may be disallowed for a given Hamil- 
tonian, corresponding to zero-valued tensor elements. 

To facilitate the description of the TRG procedure, the 
hexagonal lattice is constructed as illustrated in Fig. [2] 
At the nth step, we replace each vertex with a hexagon, 
with the initial structure denoted n = 0. We impose pe- 
riodic boundary conditions, such that the top and bot- 
tom edges are equivalent, as well as the left and right 
edges, so that the lattice effectively lies on the surface of 




FIG. 2: Construction of the hexagonal lattice used in the 
TRG procedure. Starting from the initial structure on the left 
(n = 0), at each construction step, every vertex is replaced 
by a hexagon. Periodic boundary conditions are imposed be- 
tween the top and bottom edges and between the left and 
right edges, as if the lattice is on the surface of a torus. The 
sublattices A and B are shown for the n — step. 



a torus. After the nth step, the system has N = 8-3 n sites 
and K = 4 • 3 n+1 bonds. The TRG method involves a 
renormalization-group transformation that reverses this 
construction process, mapping the system at step n to 
one at step n — 1 . 

The hexagonal lattice of any size can be decomposed 
into two sublattices A and B, such that the nearest neigh- 
bors of one type belong to the other type. As an ex- 
ample, we label the sublattices in the n — panel of 
Fig. fi] We distinguish the sublattice tensors with super- 
scripts, Ij j j or Ij fl ( j . In the partition function sum 
of Eq. ((T|), each bond index i m appears twice, once within 
an A tensor, and once within the neighboring B tensor 
linked through that bond. Thus evaluating Z consists of 
performing K tensor contractions. 

In addition to the bond variables, the general system 
we consider has quenched random degrees of freedom, 
though for notational simplicity we shall not explicitly 
show the dependence of T on these. Physical observables 
Q will be expressed as [{Q}], where (■) denotes the ther- 
modynamic average over the bond degrees of freedom and 
[•] denotes the configurational average over the quenched 
disorder. 



B. The TRG Transformation 

The TRG transformation consists of two steps, known 
as rewiring and decimation. In the rewiring step, the 
bonds of every pair of neighboring tensors T A and T B 
are reconnected, rewriting them as a contraction of two 
new tensors S A and S B . The reconnection pattern is 
illustrated in Fig. [3f a,b) and can be broken down into 
three basic cases (highlighted in different colors) involv- 
ing different orientations of the intial T A and T B ten- 
sors. In our graphical convention, the vertex where three 
solid lines meet is a T tensor and the vertex where three 
dashed lines meet is an S tensor. Indices on a tensor, 
i.e., T-4. , correspond to bonds labeled i, j, k arranged 
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FIG. 3: (Color) The TRG transformation described in Sec. Ill Bl (a) The hexagonal tensor network, with the three representative 
orientations of T A and T B tensor pairs, labeled as cases 1 through 3 and highlighted in different colors, (b) For each of the 
three cases, the rewiring step [Eq. p)]. expressing the contraction equivalently in terms of different tensors S A and S B . (c) 
After every pair of tensors is rewired, the resulting martini lattice of S A and S B tensors. The original lattice is superimposed 
in gray for reference, (d) The decimation step [Eq. l(5)]. which replaces three S A tensors by a renormalized T' A tensor (and 
analogously for S B ). (e) The final lattice of renormalized T' A and T' B tensors, with the original lattice in gray for comparison. 



counterclockwise around the tensor, with the first index 
marking the vertical bond for the T tensors and the hor- 
izontal bond for the S tensors. Thus for example the 
three rewirings shown in Fig. [3][b) denote the mathemat- 
ical identities 
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Case 1 : T A kl T B ki — ^itj^l 

k=l v=\ 
d d 2 

Case 2 : T^ m T^ = S A jiS B mi , (^) 

k=l v=l 
d d 2 

Case 3 : ^ ^Zmfc^iffc = ^fiv^miu- 

k=l v=l 

Note that the S tensors have two indices which run up to 
d (labeled by Latin letters) and one index that runs up to 
d 2 (labeled by a Greek letter). The reason why S A and 
S B must have this structure comes from the following 
derivation, which also illustrates how one can explicitly 
calculate these tensors. 

We shall describe the derivation for case 1, since 
the other two cases are analogous. The first line of 
Eq. (p]) can be expressed as a d 2 x d 2 matrix equation, 
M = S A (S B ) T , where M a0 = £ fe T^T* , S A V = Sfa, 
S B U = Sf vm . Here we use composite indices a and /3 
with d 2 states defined as a = and (3 = (m,i). As 
a real-valued matrix, M has a singular value decomposi- 
tion of the form M = UY>V T , where U, V are orthogonal 
matrices and £ is a diagonal matrix containing the d 2 
singular values of M. Once the singular value decompo- 
sition of M is calculated, the elements of S A and S B are 
given by S A V = 



where E TO is 

the vth singular value, adopting the ordering convention 
from largest to smallest with increasing v. 

After all T A and T B pairs are rewired, we have a so- 
called martini lattice of S A and S B tensors, shown in 
Fig. &c). The final step of the TRG transformation is 



decimation, which traces over the degrees of freedom in 
the triangles of the martini lattice, substituting for each 
triangle a renormalized tensor T' A or T' B . Graphically, 
Fig. [2Jd) shows the decimation of three S A tensors to 
form T' A and of three S B tensors to form T' B . The 
corresponding expressions in terms of tensor components 
are 



EqA qA qA rplA 
D yjl D l-ih D hjS — 1 1/7(5' 

j,l,h=l 
d 

EQB qB qB _ rpll 

^ i/mi i'yh hmS ^7 

m,i,h— 1 



(3) 



-ilB 



The final renormalized tensor network of T' A and T' B is 
shown in Fig. G^e). 

The partition function Z , a contraction over all bonds 
connecting the tensors, Eq. ([1]), is exactly preserved 
through this transformation, as the hexagonal lattice is 
coarse-grained from a step n to a step n — 1 structure. 
However, the indices of the renormalized tensors run from 
1 to d 2 instead of 1 to d, so that if the TRG were iter- 
ated, arbitrarily large tensors would result, making nu- 
merical implementation difficult. This problem is related 
to a general feature of position-space renormalization on 
lattices: except for specially tailored geometries (i.e., hi- 
erarchical lattices [18|-|20j), the number of couplings in 
the renormalized Hamiltonian grows with each coarse- 
graining. For the TRG, we can tackle this issue in a sys- 
tematic fashion by truncating the index range with an up- 
per bound D. In Eq. © for T' A and T' B ,_ we shall allow 
the indices v,"f, and 5 to run only up to d = min(d 2 , D). 
This is equivalent to using truncated matrices S A and 
S B in the rewiring step, where S is the first d columns 
of the d 2 x d 2 matrix S A and S B is the first d columns 
of S B . As a result, the rewiring becomes approximate, 
M w S A (S B ) T . But since the first d columns correspond 
to the largest singular values, the approximation is rel- 
atively accurate even for small D and rapidly converges 
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FIG. 4: (Color online) Duality mapping between spin states 
on the triangular lattice and bond variables in the tensor net- 
work. The variables Si = ±1 at the triangle corners corre- 
spond to Ising spins in the Hamiltonian of Eq. (Q. The bond 
variables Oi are products of the Sj connected by the bond. Up 
and down triangles yield type A and B tensors respectively. 

as D is increased 0, Q . With this cutoff, the maximum 
size of the tensors is bounded as the TRG procedure is 
iterated and we can extract numerically thermodynamic 
information from flows within a finite-dimensional space 
of real- valued tensor elements. 



III. TRG FOR QUENCHED RANDOMNESS: 
THE BOND-DILUTED ISING MODEL 

A. The Bond-Diluted Ising Hamiltonian and Its 
Mapping onto a Tensor Network 



a T B , as shown in Fig. |4j For spin variables s,, Sj, s& 
in a given triangle in the manner illustrated in the fig- 
ure, we define corresponding edge variables a m as the 
products of neighboring s TO , i.e., for the type A triangle, 
ci = S3S1, a 2 = s%S2, 03 = S2S3 and for the type B tri- 
angle, cti = sis 2 , cr 2 = s 2 s 3 , <7 3 = S3S1. Since s m = ±1 
and a m = ±1, we can now introduce a composite index 
i m = (5 — <J m — 2s m )/2 which runs from 1 to 4 and de- 
scribes the four possible states of the rath triangle edge. 
Letting J m be the coupling J i3 associated with this edge 
and H m = H^ be the edge magnetic field, then the ten- 
sors for the two triangles types are: 

TLi 3 =e*^ ^-+^( 1 +-)-)p(a 1 a 2( x 3 ) 

• P(criSlS 3 )P((T2S2Sl)P((T3S3S2) , 

(6) 

■ P{criSis 2 )P((T2S2S S ,)P(a3S 3 Si) , 

where P(x) = (1 + x)/2 is a projection operator. The P 
factors in the tensors remove the bond states that do not 
correspond to a physically allowable spin configuration. 
As a result of the projection operators, only 8 out of 
the 64 elements in the tensor are nonzero. These are 
listed, for the first rcnormalization step, in the 3rd and 
6th columns of Table |l] for T A and T B respectively. 



The general Hamiltonian for a quenched random Ising 
system is 

-pU = Y^ [J l3 s tSj + Hijisi + Sj )] , Si = ±1 , (4) 

where /3 = 1/ksT, Jij and Hij are respectively the lo- 
cal spin-spin coupling and magnetic field for sites i and 
j, and (ij) denotes a sum over nearest- neighbor pairs of 
sites. Although this Hamiltonian encompasses a va riety 
of models, all the way to the random- field spin glass [2l| . 
we shall here focus on a the bond-diluted Ising case, 
where the interaction constants are distributed with 
a quenched probability V(Jij) of the form 

V(Jij) = P S(J l3 ) + p)8{Jij - J). (5) 

Here J > 0, implying ferromagnetism, and p is the frac- 
tion of missing bonds. While we restrict our attention to 
the zero magnetic field subspace, = 0, formally the 
local fields will be kept in the Hamiltonian in order to 
take derivatives to obtain thermodynamic functions. 

Starting with the Hamiltonian of Eq. ((4]) on a trian- 
gular lattice, a duality transformation allows us to ex- 
press the partition function as a hexagonal-lattice ten- 
sor network. (The duality for Potts spins would gen- 
erate three-point interactions, which would be included 
in the definition of the tensor Ti 1 i 2 i 3 .) Each triangle in 
the triangular lattice corresponds to a tensor, with up 
triangles associated with a T A and down triangles with 



Local Magnetization and Spin-Spin Correlation 
Function 



In order to derive expressions for thermodynamic 
quantities in the tensor formalism, let us now restrict 
the notation T A and T B to tensors in the zero magnetic 
field subspace. We place a local magnetic field Hk only 
at a single location k. Let us call the two tensors which 
share this bond T A and T B . These are the only two ten- 
sors in the system whose components are modified by the 
local field. The corresponding partition function is 

J7 \ * rpA rpB i _ | rpA J^B _ _ _ rpB 

ii,...,i K 

(7) 

Without loss of generality we take the contraction of the 
f A and f B tensors to be Case 2 in Eq. (0), since the 
derivation proceeds analogously for the other Cases. 

The local magnetization is ra^ = ((s, + Sj)/2) = (Sk) 
for the sites i, j associated with the bond k. In terms 
of the local magnetic field -ff/c, the magnetization ra^ is 
given by the derivative 



1 dlnZ 



H k =Q 



2 dH k 

— Y \T A 



12*3 i*^5*3 ' ' ' *kilim Wnio 



. (8) 



I rj^A rpB 



T I)" 
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11*2*3 


ttt 


in 
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■M- 


421 


e i(-J!-J 2 + J3+2H3) 







412 


e i(-J!+J 2 -J3 + 2H 2 ) 





m 
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ut 
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e 5(Jl-J2-J3-2H!) 


_o 2 (Jl-J2 — J 3 ) 




333 


e i(Jl + J 2 + J 3 -2H 1 -2H 2 -2H3) 




+- J 3 ) 
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e ±(Jl+J 2 + J 3 -2H 1 -2H 2 ~2H 3 ) 
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TABLE I: The tensor elements for the bond-diluted Ising model, as defined in Sees. IHIAl and lTlI Bl for the first renormalization 
step. The first column gives the spin state (si,S2,S3) for a triangle of the original triangular lattice, following the convention 
of Fig. [J For the type A triangle, the next three columns show the associated composite indices (11,12,13), and the tensor 
elements T A i2i3 , D A i2i3 . The last three columns show the analogous information for the type B triangle. All tensor elements 
not shown are zero. 



where the differentiated tensors are 



C. Details of the Numerical Implementation 



D] 



dT A - ■ 

*fc*!*r, 



dHi, 



D 



dT B 



Hi,=0 



Hi k =0 

(9) 



The nonzero elements of D A and D B are shown, for the 
first renormalization step, in the 4th and 7th columns of 
Table H 

After taking the average over the disorder, the first 
and second terms in the brackets on the right-hand side 
of Eq. (|5|) are equivalent, so that 

[mic] =[(Sk)] = 

y-\ rpA rpB r)A rpB 

/ j *i* 2 *3 14*5*3 ' " ' ikilim ikinio 
ii,—,*K 

(10) 

A similar derivation for the correlation function yields 



[(SkSi)} = 



ry — 1 \ rjiA rpB 7~\ j4 ryiB 



■ • • D A - T B ■ 



(11) 



We shall be interested in long-range correlations, as an 
indicator of thermodynamic phase behavior. In this case, 
the four individual Si spin-spin correlations that make up 
the [(S/ciS;)] are approximately equal: [(SkSi)} 



where Si is either of the spins contributing to Sk and Sj 
is either of the spins contributing to Si. Hence we shall 
use [(iSfe.!?;)] and [(siSj)] interchangeably in the rest of the 
text. 



To calculate the long-range spin-spin correlation func- 
tion [(S'fcS'i)], we start with a finite hexagonal lattice after 
n construction steps, with size varying between n = 7—10 
steps (N = 17496 - 472392 tensors). The bonds k and 
I are chosen to be at the maximum separation within 
the lattice, taking periodic boundary conditions into ac- 
count. For a given realization of the disorder, the sum 
on the right-hand side of Eq. (fTT|) is evaluated by doing 
n TRG transformations, which yields the contraction in 
terms of four renormalized tensors in the n = structure. 
These last four tensors are directly contracted. A similar 
process yields the value of the partition function Z which 
is the denominator in Eq. ([lip. The configurational av- 
erage is taken over 200 — 300 realizations, implemented 
by randomly assigning the on the initial lattice ac- 
cording to the probability distribution in Eq. ([3]). The 
tensors on the original lattice, i.e., in Eqs. (jSJ) and ©, 
have index range d = 4. For subsequent tensors, we use 
a cutoff parameter D = 8 — 14. 



Some tensor elements tend to grow exponentially in 
magnitude as the TRG transformation is iterated, which 
poses potential numerical difficulties. To counteract this, 
we take advantage of the fact that we can always fac- 
tor out a constant from each tensor without changing 
the physics. For each tensor during each TRG iteration, 
the factor extracted is equal to min(T max , 2) where T max 
is the maximum absolute value of the tensor elements. 
Keeping an upper bound of 2 on this extracted factor 
slows down the decay of most tensor elements to zero, 
which would otherwise lead to other numerical artifacts. 
We keep track of the total extracted factors in the numer- 
ator and denominator of Eq. pip , which are then used 
in calculating the final correlation function value. 
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FIG. 5: (Color online) The long-distance spin-spin correlation 
[(sjSj)] as a function of temperature 1/J, calculated using the 
TRG method for bond dilution probability p = 0.1 and cutoff 
parameter D = 8. The curves for four different initial tensor 
network sizes N are shown. 
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FIG. 6: (Color online) The long-distance spin-spin correlation 
[(s»Sj)] as a function of temperature 1/J, calculated using the 
TRG method for bond dilution probability p — 0.1 and net- 
work size N = 157464 tensors. The curves for four different 
cutoff parameters D are shown. 



N using cutoff D = 8, while the latter varies the cut- 
off D at fixed size N = 157464. Away from the critical 
temperature, where widely separated spins are uncorre- 
lated, [(siSj)] « [(s^) 2 ], and wc expect distinct limiting 
behaviors for the two different thermodynamic phases in 
the system: at low 1/J in the ferromagnetically ordered 
phase [(siSj)] — > 1, while at high 1/ J in the paramagnetic 



phase [(SiSj 



0. The temperature region where one 



sees a smooth transition between these two regimes for 
finite systems, illustrated in Figs. [5] and HI gives a rough 
indication of the phase transition temperature 1/ J c in the 
thermodynamic limit. With increasing N in Fig. [5] and 
increasing D in Fig. [6j the transition becomes sharper, 
as our truncations converge toward the exact result for 
an infinite system. The probability p = 0.1 at which 
these results are calculated is smaller than the thresh- 
old value p c w 0.653 [22| . above which the triangular 
lattice no longer percolates. For p > p c we would not 
see a transition region: the paramagnetic phase exists 
at all temperatures, since islands of ordered spins of size 
~ O(N) become exponentially improbable. 

To obtain an accurate estimate of the exact transition 
temperature 1/J C , we can employ the following finite- 
size scaling relation, which describes the ratios of the 
correlation functions at three different system sizes N±, 
N 2 , and N 3 when J = J c 



In 



g(N 3 ) 

g(N 2 ) 



In 



(f) 



(12) 



where g(N) is the long-distance correlation function 
[(sjSj)] for network size N. For the ith system, at the 
temperature region where g(Ni) decays rapidly to zero 
(J just smaller than J c ), the decay is approximately ex- 
ponential in J, 



ln(5(7Vj)) w AiJ - Bi 



(13) 



for some constants Aj and Bi. This exponential behav- 
ior for three different system sizes is shown in Fig. [7] for 
p = 0.25 and 0.55. To calculate Ai and Bi, we do a 
weighted linear least squares fit to \n(g(Ni)) vs. J data 
in a region of J where the relative uncertainty (from the 
configurational average) for the data points is less than 
15%. Plugging Eq. (13]) into Eq. (12]) with J = J c , we 
can solve for J c in terms of the Aj, Bi, and Ni, 



Jc 



(B 2 -B 1 )bx(^ 



(S 2 -B 3 )ln(f) 



(A 1 -A 2 )\n(^)+(A 3 



(14) 



IV. RESULTS 



Representative results for the long-distance spin-spin 
correlation function [{-SiSj}] as a function of temperature 
1/J at bond dilution p — 0.1 are given in Figs. [5] and H] 
The former shows curves for various tensor network sizes 



Carrying out this calculation across the entire p range 
for Ni = 17496, N 2 = 52488, and N 3 = 157464 at 
D = 12, we obtain the phase diagram shown in Fig. [T] 
For comparison we also plot the same phase diagram 
obtained from a rigorous approximation scheme for the 
bond-diluted Ising model free energy [l7| . which can 
be considered exact on the scale of the figure. The 



7 



c ■ 
o 

4-» _ 

u 

£Z 

3 . 

c 
o 

■4-» 

(U 

L_ 
i_ 

o 
u 

01 
cn 

fD 
I 

o> 

£Z 

o 




-5 
-10 
-15 
-20 
-25 

-1 

-2 
-3 
-4 
-5 
-6 




0.370 



0.375 



0.380 



0.385 









□ 

□ 




• N = 17496 
t N = 52488 
□ N = 157464 



0.80 0.82 0.84 0.86 0.88 
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FIG. 7: (Color online) Data points show the log of the long- 
distance spin-spin correlation, ln[(siS 3 -)], as a function of inter- 
action strength J for three different system sizes N and two 
different bond dilution probabilties p (top panel: p = 0.25, 
bottom panel: p — 0.55). The weighted least squares lin- 
ear fits, shown as solid lines, yield the coefficients At and Bi 
in Eq. (|13[) . which allow one to estimate J c through finite 
size scaling, Eq. <|14|) . The resulting values of J c are 0.388 
(p = 0.25) and 0.871 (p = 0.55). 



CONCLUSIONS 



We have shown how the TRG approach provides an ef- 
ficient and precise method for calculating thermodynamic 
properties of a quenched random classical model — the 
triangular-lattice bond-diluted Ising Hamiltonian. By 
expressing the partition function and related quantities 
such as spin-spin correlation functions in terms of tensor 
networks, they can be readily evaluated through TRG for 
large lattice sizes. In combination with finite-size scaling 
ideas, the result is a precise estimate of the phase dia- 
gram. If desired, convergence to the exact critical prop- 
erties can be achieved by increasing the cutoff parameter 
defining the index range of the tensors. 

The bond-diluted Ising model is only a first step in 
the exploration of disordered systems using TRG: the 
methods presented here are easily extended to frustrated 
Hamiltonians exhibiting spin-glass behavior and the re- 
sulting complex multicritical phase structures. The nu- 
merical accuracy of the technique will be a valuable fea- 
ture in probing analytical conjectures on the exact loca- 
tions of spin-glass multicritical points [25l - [28| . 
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